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Abstract 


In order to tackle parameter estimation of photocounting distributions, polykays of acting 
intensities are proposed as a new tool for computing photon statistics. As unbiased estimators 
of cumulants, polykays are computationally feasible thanks to a symbolic method recently de¬ 
veloped in dealing with sequences of moments. This method includes the so-called method of 
moments for random matrices and results to be particularly suited to deal with convolutions 
or random summations of random vectors. The overall photocounting effect on a deterministic 
number of pixels is introduced. A random number of pixels is also considered. The role played 
by spectral statistics of random matrices is highlighted in approximating the overall photo¬ 
counting distribution when acting intensities are modeled by a non-central Wishart random 
matrix. Generalized complete Bell polynomials are used in order to compute joint moments 
and joint cumulants of multivariate photocounters. Multivariate polykays can be successfully 
employed in order to approximate the multivariate Poisson-Mandel transform. Open problems 
are addressed at the end of the paper. 

Keywords: photocounter, mixed Poisson distribution, non-central Wishart random matrix, 
symbolic method of moments, cumulant, polykay. Bell polynomial 

1 Introduction 

The recent renewed interest in multivariate photocounting originates within astronomical litera¬ 
ture in connection with extrasolar planet detection methods |l3] or, more in general, in speakle 
patterns produced by direct imaging m- A photovent occurs when light striking a pixel causes 
one or more electrons to be ejected. Within fixed time intervals T, multivariate photocounters 
are modeled by non-negative random vectors counting photoevents from a set of d pixels. Their 
stochastic fluctuations depend on intensities of light, encoded in a vector I = (Ii,... ,Id), with 
joint distribution /r(d/) on The photocounters {Ni,i = 1,... ,d} are assumed condition¬ 

ally independent and distributed according to a Poisson law parametrized by I. The multivariate 
Poisson-Mandel transform m 



( 1 . 1 ) 


gives the joint distribution of AT = (W,..., Nd) when k = (fei,..., kd) is a multi-index of non¬ 
negative integers. The availability of closed form formulae for (ffAD depends upon the stochastic 
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model fi. For a wide range of univariate stochastic intensities, these formulae have tractable ex¬ 
pressions and sufficient statistics can be recovered by means of likelihood methods [3]. In [29], an 
overview is given on the methods available to analyze the data sampled in photocounting, as for 
example the photon counting histogram (PCH). Fluorescence cumulant analysis (FCA) is indi¬ 
cated as the first theory that describes the effect of sampling time by measuring the spontaneous 
intensity fluctuations of fluorescent molecules |20|. This measurement is done by using facto¬ 
rial cumulants. Let us recall that if M(t) denotes the moment generating function of a random 
variable (r.v.), then cumulants are the coefficients of log(M(t)) and factorial cumulants are the 
coefficients of log(M(log(l -|- t))). Cumulants have special properties and if some of these proper¬ 
ties can be deduced from data, then special stochastic models can be inferred m- For example, 
if conditioned cumulants linearize, then the recorded data can be modeled by using (ll.lll . The 
same relationship holds for factorial moments. 

Inspired by FCA, the main goal of this paper is to extend cumulant analysis to the multivariate 
case, which is still a challenging problem in particular when I are entries of a random matrix 
m- The preliminary contributions given in the literature indicate factorial moments as simpler 
expressions to be used mEi]. 

In working with random matrices, the method of moments is still extensively used [Ij. In 
general, some conditions on moments or moment generating functions need to be required and 
applications of this approach to random matrices rely on some advanced combinatorial tools as 
for example zonal polynomials or hyper geometric functions [2a|25l[27j. A first way to overcome 
this drawback is free probability, a non-commutative theory of probability |2T]. Indeed, random 
matrices are non-commutative objects whose large-dimension asymptotic are usually analyzed by 
using free probability. However, there are some aspects of random matrix theory to which the 
tools of free probability are not sufficient by themselves to resolve. 

In this paper, we propose the employment of a method, called the symbolic method of moments 
[8], which can be considered the commutative counterpart of free probability. In the symbolic 
method of moments, a sequence of numbers is represented by a symbol, called umbra, through a 
linear operator, sharing many properties of expectation. The elements of the sequence are called 
moments. This symbolic approach overcomes the well-known moment problem, since a sequence 
of numbers is dealt as it was a sequence of moments with no reference to any probability space 
as within free probability. It is mainly a tool to perform computations: the matching with r.v.’s 
is done a-posteriori. In difference from free probability, just one operator is employed within the 
symbolic method of moments but the same sequence of moments may correspond to more than 
one symbol. For a pair of commutative r.v.’s, freenes£| is equivalent to claim that at least one 
of them has vanishing variance. So freeness is a pure not-commutative device that is why the 
symbolic method of moments may be considered an analogous of free probability in a commutative 
field. One of the strengths of this method is that questions on convergence of moment generating 
functions may be discarded and only polynomials involving umbrae are employed. 

That having been said, the novelty of this paper is twofold. The first has implications beyond 
photocounting and involves the employment of the symbolic method of moments within random 
matrices. The second novelty is strictly related to photocounters and their statistics. In FCA, 
factorial cumulants are employed and calculated from the moments of the recorded photon counts. 
Quite recently efficient algorithms [9] have been developed in order to compute experimental 
measurements of cumulants, known in the literature as polykays. The sampling behavior of 
polykays is much simpler than sample moments m but their employment was not so widespread in 

^In free probability, free r.v.’s correspond to classical independent r.v.’s. 
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the statistical community, due to the past computational complexity in recovering their expression. 
When complex amplitudes of incoherent waves have independent circular Gaussian distribution, 
the resulting acting intensity I is the diagonal of a non-central Wishart random matrix and we 
show how spectral polykays [TO] can be fruitfully employed in order to estimate their cumulants 
and then factorial cumulants of photocounting. 

We prove that photocounter cumulants have a simpler expression compared with moments and 
factorial moments, taking advantage from the plainness of cumulants of the non-central Wishart 
random matrix. Thanks to the symbolic method of moments, the generalization to multivariate 
framework is straightforward. And this is an additional novelty of the paper since far fewer results 
can be found in the literature on compound Poisson random vectors [22|, due to the difficulty of 
managing their distributions. 

Since the symbolic method of moments is a new tool for photocounting, we get the opportunity 
of introducing this theory by dealing with a new measure describing the overall photocounting 
effect, together with its generalization to superposition of a random number of incoherent waves. 
New formulae are proposed to perform all these computations. Implementations of these formulae 
in Maple are available on demand. 

The paper is organized as follows: the symbolic method of moments is sketched in Section 
3 for the univariate case and in Section 4 for the multivariate case. Section 2 shows how to 
model acting intensities when complex amplitudes of incoherent waves have independent circular 
Gaussian distribution. Polykays and their properties are recalled in Section 3. Formulae for 
computing moments, factorial moments and factorial cumulants of multivariate photocounters are 
given in Section 4. A series expansion of the multivariate Poisson-Mandel transform shows how 
to by-pass the multi-dimensional integral in (jl.l|i and to use multivariate polykays for numerical 
approximations. Some open problems are suggested at the end of the paper. 

2 Acting intensities modeled by Wishart random matrices 

The complex amplitude ijj{x,y) of a wave can be modeled as ijj{x,y) = m(x,y) + X{x,y) [I], 
where m{x, y) is a deterministic term proportional to the wave amplitude without turbulence, and 
X{x,y) is a random term distributed according to a zero mean complex Gaussian distribution. 
This choice depends on the central limit theorem, since X{x, y) represents the uncorrected part of 
the wave amplitude caused by errors. The instantaneous intensity is defined as I{x, y) = \X{x^ y) + 
m(x,y)p. When p incoherent waves are considered, their superposition field intensity /[p](x,y) is 
obtained by summing the intensities of each wave, that is I[p]{x,y) = , y)\‘^ ■ 

A multidimensional framework is necessary when d pixels are involved and correlations among 
wave amplitudes at d different positions are encoded in a full rank covariance matrix S. Then the 
complex amplitude of the i-th incoherent wave at (x, y) is given by y) = mi{x, y) + Xi{x, y), 
where the deterministic term mi[x,y) is a d-dimensional vector and the random term Xi[x,y) 
is a d-variate circular complex Gaussian random vector with zero mean and full rank Hermitian 
covariance matrix S. By omitting the notation (x, y) for brevity, the resulting vector of intensities 
hv] = {h,[p]i ■ ■ ■! ^d,[v]) has components 

p 

^j,lp] = '^\^ij for j = 1,2,... ,d (2.1) 

i=l 
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with Xij = {Xi)j and rriij = {mi)j. The elements of are on the diagonal of a non-central 
Wishart square random matrix of order d 

p p 

“1“ With ]\d — ^ ^ TTl^TTLi^ ^2.2^ 

i=l i=l 

usually denoted by Wd{p)- In equation (12.21) . f denotes the conjugate transpose. In the literature 
Q = Tj~^M is called the non-centrality matrix. 

The intensity /[pj is the parameter of a photocounting vector = (A^i,[p], • • • ,IVd,[p]) with 
Poisson distribution (jl.ip . The Poisson r.v.’s IVijpj,... ,IVd,[p] are conditionally independent, that 
is P (ATW = k I 4j) = n •=! P = kj I Ijm) ■ 

By recalling that complex Gaussian random vectors of dimension d can be replaced by real 
Gaussian random vectors of dimension 2d, whose hrst d-components represent the real part and 
the last d-components represent the imaginary part, an alternative way to compute Ij,[p] in (12T]) is 
^j,lp] ~ [lIe(Xjj -|- rriijY + Ini(Aijj -|- , for j = 1, 2,..., d. However, this formula picks 

out elements on the diagonal of a real non-central Wishart random matrix of dimension 2 d and 
therefore is less efficient. Then, we refer to complex non-central Wishart random matrices (|2.2p . 


3 Overall photo counters 

Let us introduce the notion of overall photocounter. 

Definition 3.1. The overall photocounter is A/jp] = IVi,[p] + • • • + IV^jp] if P incoherent waves hit 
d pixels and A^qp] denotes the number of photoevents of the j-th pixel labeled with j = 1,... ,d. 

Since the convolution of two (or more) mixed Poisson distributions is itself a mixed Poisson 
distribution, with mixing densities the convolution of the two (or more) mixed Poisson distri¬ 
butions m, then the overall photocounter A/'q] has a mixed Poisson distribution with random 
parameter /i,[pj H-h IdM = Tr [Wdip)] ■ 

Proposition 3.2. If (•)j denotes the lower factorial and S{i,k) are Stirling numbers of second 
kind, then 


A[(A/-[pj)J =A{(Tr[IT,(p)])‘} 


and 


E 


(A/'w)* 


^S{z,k)E[{Tr[Wd{p)])’^}. 

k=l 


Proof. Conditioned factorial moments of Poisson r.v.’s are powers of the random parameter, that 
is E [(A/'[p])J di,[p], • ■ ■, Id,\p]] = (Tr [Wd{p)]y ■ Then factorial moments follow by taking the overall 
expectation. Moments follow by taking the overall expectation of x® = '^’(b ^)(2;)fc) after 

having replaced the indeterminate x with A/Jq. □ 

The distribution of the overall photocounter can be computed as 


p (A/'w = ^) = ^ E ^ {(Tr [Wdip)]?^'} = ^ E ^pn+.(ci,..., (3.1) 

2=0 2 = 0 

with Yk{xi,... ,Xk) (complete) exponential Bell polynomials [12], see also Section 6.2, and {ck} 
cumulants of Ti \Wd{p)\. In [6|, moments of Th:\Wd{p)] have been expressed by using integer 
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partitions d and cyclic polynomials d- Let us recall their expression to highlight their complexity 
with respect to cumulants. If Cfc(S) denotes the /c-th cyclic polynomial Cfc(si,..., Sk) with Sk = 
Tr[S^], then 

k 

E[T,[Wd{pt] = k\Y, 

i=o 

where TrA(MS) = njg{i 2 }[Tr(MS'^“^)]''^ and Ca(S) = njg{i 2 algorithm to 

compute (13.2p is available in [6j relied on the symbolic method of moments. 

The superposition of incoherent waves is fully employed when Tr [ITd(p)] is written as 

dp p 

Li,[p] H h Id,{p] = EE \Xij + mjjp — \Xii + TTiiip + • • • + \Xici + midp, (3-3) 

j=l i=l i=l 

since (|Xii+miip,..., l^id+midp) are independent row vectors for i = 1,... ,p. In difference from 
cumulants, neither moments nor factorial moments of take advantage of the decomposition 
on the right hand side of equation (13.31) . Instead, conditioned cumulants linearize and are equal 
to the random parameter of the overall photocounter for all non-negative integers k, that is 

d 

Cumfc(iVi,[p]-k--- + lVd,[p]|li,[p],...,/d,[p]) = ^Cumfc (iVj,[p]|4[p]) = Tr[ITd(p)]. (3.4) 

i=i 

Unconditioned cumulants will be computed in the next section, since the symbolic method helps 
in shortening the proofs. Here, we limit ourselves to observe that 

p 

Cumfc (A/jpj) = ^ Cumfc [iVi* H-h Ndi] (3.5) 

i=l 

with Nji {i = 1,... ,p, j = 1,... ,d) the photocounter related to the j-th pixel and the i-th wave. 
Approximations of the distribution in ()3.1I1 involve cumulants of Wishart random matrices which 
have a plainer expression compared with moments (13.2p . as we will show in the next section. 




ai-7 


m(A)! 


E 

Xhk-j 


p 




m{xy. 


Ca(E) 


, 


(3.2) 


4 The symbolic method of moments 

In the symbolic method of moments, an alphabet A = {a,/3,7,...} of indeterminates, named 
umbrae, is considered and any umbra is related to a complex number sequence {a^} by a suitable 
linear functional E. The functional E : C[^] —>■ C is defined on the polynomial ring C[Al], and 
such that E[a*^] = for all non-negative integers k > 1. We assume E[l] = 1 so that uq = 1. 
The sequence {a^} is the sequence of moments of a and we say that {a^} is umbrally represented 

^Recall that a partition of an integer fe is a sequence A = (Ai, A 2 ,..., At), where \j are weakly decreasing 
integers and Aj = fe. The integers Aj are named parts of A. The length of A is the number of its parts and 

will be denoted by /(A). A different notation is A = (1’'^, ...), where rj is the number of parts of A equal to 

j and ri -I- r 2 -I- • • • = ((A). We use the classical notation A h fe with the meaning “A is a partition of fe”. Set 
m(A) = (ri, r 2 ,...) and m(A)! = ri!r 2 ! • • • . 

^The i-th cyclic polynomial is Ci{x \,..., Xi) = ■ ■ ■ x^ with d'x = i!/(U'^ri!2’"^r2! • • • )• S®® [S] for more 

details. 
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by a. Two umbrae can represent the same sequence of moments, that is E[a^] = E[ 7 ^] for all 
non-negative integers A; > 1. In such a case we said that a is similar to 7 , in symbols a = 7 . The 
operator E factorizes on distinct umbrae, that is E[a*/3-^ • • • 7 ^] = E[a®]E[/3'^] • • • E[ 7 ^] (uncorrela¬ 
tion property). A conditional evaluation has been introduced in [11] satisfying E[q!*/3'^ • • • 7 ^|a] = 
a*E[/3^']---E[7*^]. 

Special umbrae are: 

a) the unity umbra u whose moments are { 1 }; 

b) the augumentation umbra e whose moments are E[e^] = (5o,fc) for all non-negative integers k 

with (5 o,A: the Kronecker Delta; 

c ) the Bell umbra whose moments are the Bell numbers |9] ; 

d) the singleton umbra whose moments are E[y;^] = for all non-negative integers k. 

The sequence of cumulants {ck} of a is defined a^ 

k>l \ k>l J 

Definition 4.1. If A is an integer partition and {a^} (respectively {cfc}) is the sequence of 
moments (respectively cumulants) of a, the product a\ = • • • (respectively c\ = C 1 C 2 " 

is said associated to the partition A. 

For example, if A = (1^,3) then ax = 0 ^ 03 . Auxiliary umbrae are introduced as special 
symbols representing operations among moments. For example, summations of n distinct but 
similar umbrae have moments 

u\ 

E[(a + • • • + a'f] = g dx CA, where dx = ( 2 !)r- 2 ^^!... (4-2) 

and cx is the sequence of cumulants of a associated to the partition A. We denote by n.a the 
auxiliary umbra representing the sequence of moments (14.2p . A generalization of n.a is the 
auxiliary umbra 7 . 0 ;, obtained from n.a by replacing n with 7 . If {< 7 ^} is umbrally represented by 
7 , this replacement corresponds to replace {n^} with {gk} in (14.2h . that is 

E[(7-a)^] = X] S'hA) CA- (4.3) 

Ahfc 

Since the dot corresponds to a summation, 7.0 represents a symbolic summation 7 times of the 
umbra a. The auxiliary umbra 7.0 is named the dot product of 7 and a and is the symbolic 
counterpart of what we call generalized random sum. 

Definition 4.2. The r.v. with sequence of moments (|4.3p is named generalized random sum. 

Suitable choices of 7 and a correspond to suitable choices of {gk} and {c^} in ()4.3p and allow us 
to recover moments of special (auxiliary) umbrae. For example jd.y = X*/3 = u- Two special dot- 
products need to be mentioned separately: the a-cumulant umbra x-ct, representing the sequence 

^Within formal power series, equation gu holds independently from questions of convergence |26| . 
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of cumulants {cfc} in (|4.1I) . and the a-factorial umbra a.x, representing the sequence {fk} such 
that E[(Q;)fc] = fk for all non-negative integers k. By analogy with r.v.’s, the complex numbers 
{fk} are said factorial moments of a. Therefore the umbra (x.a).x = x-(<a-x) represents the 
sequence of factorial cumulants. The following definition states when an umbra may be replaced 
by a r.v. 

Definition 4.3. An umbra a represents a r.v. X, if a umbrally represents the sequence of 
moments {E[X^]}, that is E[a^] = E[X^] for all non-negative integers k. 

In particular, the dot-product n.a represents the summation of n i.i.d.r.v.’s. The auxiliary 
umbra 7.0 represents a generalized random sum. If moments are defined only up to some non¬ 
negative integer k, then sequences of only k elements are considered. 

Symbolic representation of overall photocounters. Assume to denote by co^ the Wishart umbra 
representing the outer product of Gaussian vectors (X* -|- Tnj)l(Xj -|- m*). From equation (12.21) . 
Tr \Wd{p)] is a summation of traces of (Xj -|- mj)^(Xj -|- m*). Therefore the umbra ijJd.yp] = -|- 

• • • represents Tr \Wd{p)] ■, with the subscript i = 1,... corresponding to the subscript of 

the mean vector m^. In [ 6 ], a different symbolic representation of Tr [IFd(p)] has been provided in 
order to speed up the implementation of formula (|3.2I) and to take advantage of the non-centrality 
matrix. This symbolic representation does not take into account the additivity property of traces, 
which instead is very useful in dealing with mixed Poisson distributions. Since a Poisson r.v. with 
random parameter A is represented by the umbra 7 ./?, with 7 the umbra representing moments 
of A [12] then the overall photocounter is represented by 

+ • • • + (4.4) 

where the right hand side of (14.4p is obtained from the left distributive property of the summation 
with respect to the dot-product m- If rrii = m for i = 1,... then uJd,\p\ = P-oJd- 

Proposition 4.4. The umbra represents the sequence of factorial moments o/A/jp]. 

Proof. The sequence of factorial moments of A)p] is represented by ooa,[p]-/3.x- The result follows 
by observing that 13.x = u and LOa,[p].u = □ 

Thanks to the symbolic representation of A/’jp], cumulants of A/jp] can be computed by using 
cumulants of Tr [Wdip]] . 

Theorem 4.5. Cum^ (A/jp,) = S{k,i) [p{i - l)!Tr(E0 - dTr(ME*-i)] . 

Proof. Cumulants of A/jp] are represented by x-(wd,[p]-/3) = (X'Wd,[p])-/3- The result follows from 
dO]), by observing that E I [(x.Wd,[p])./3]''I = Yli=iS{k,i)E [(x-w^jp])*] (cf. [H]) andE [(x-w^.m)*] 
is the z-th cumulant of Tr [ITd(p)] • Its expression is given in [^. □ 

From Theorem 14.51 the additivity property (|3.5p of cumulants can be recovered since 

Cumfc(Tr[IFrf(l)]) = (A;-1)!TV(S'') - A:!Tr(m|miE*’-^) i = l,...,p. (4.5) 

Factorial cumulants FCum^ (A/Jpj) of A/)p] are equal to cumulants of Tr [IF(i(p)] . 

Theorem 4.6. FCumfc (Xjp]) =p{k- l)!Tr(S^) - k\Tr{ME’^-^). 

Proof. Factorial cumulants of A/jpj are represented by x-{'^d,[p]'(3).x = X-^d,ip] since /3.x = u- 
Moments of x-Wd,[p] are cumulants of Tr [IFd(p)] . □ 
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4.1 Randomized overall photocounting effect 


In literature on photocounting effect, the number of incoherent waves hitting the pixels has been 
considered deterministic. Here we assume this number described by a r.v. P. 

Definition 4.7. The randomized overall photocounting effect A/jpj = + • • • + is the 

number of multivariate photoevents, if P incoherent waves are superimposed on d pixels with 
intensity 

p 

Ij,P = '^\Xij + for j = 1,2,... ,d. (4.6) 

i=l 

Note that equation (14.6p is obtained from (12.11) . setting rriij = rrii for all non-negative integers 
i, in order to have Xi + m ~ N{m, S) and m = (mi,..., m^). From (j4.6p . the random parameter 
of iVi,pH- \-Nd^p may be written as/i,pH -|-/d,p = Tr[(Xi-|-m)i(Xi-|-m)] = Tv[Wd{P)]. 

Proposition 4.8. If Af^p] is the randomized overall photocounter and p is the umbra representing 
the r.v. P, then 


E 




= E 


{p. UJd.P)’' 


E[{Afi 


IP] Ik] 


= E 




Cumfc (A/^pj) = E {{x.p).{uJd-P)}’' , FCumfc (A/jp]) = E {x-P-uJd) 


Proof. From Dehnition 14.71 the randomized overall photocounter Afip] is obtained from A/jp] by 
replacing p with P. From (14.41) the umbral counterparts of A/jpj and A/jpj are p.uJd-P and p.cvd-ld 
respectively, since ujd,[p] = P-^d- Factorial moments are represented by p.ujd-l3-x and the result 
follows as 13.x = u. Cumulants are represented by x-P-^d-P and the result follows from the 
associativity property. Factorial cumulants are represented by x-P-^d-l3.X = X-P-^d- D 

Corollary 4.9. //A/jpj is the randomized overall photocounting effect, then J\f[p] = Xlili-^*,[ 1 ]) 
where represents the overall photocounter of the i-th wave. 

Proof. Since p.{ujd.l3) = {p.L0d).l3, the result follows observing that these two umbrae represent 
respectively -^IP]- equality in distribution follows since they both have con¬ 
vergent moment generating function. □ 


In the following, when no misunderstandings occur, we denote simply by A^[i]. By using 
the symbolic method, moments, factorial moments and cumulants of the randomized overall 
photocounter A/’fpj can be easily recovered. 


Proposition 4.10. If Afip] is the randomized overall photocounter, then 


i) E 




Ex^kdxE[P^W] Cum,(iV[,); 


n) E [(A/-[p,),] = ExrkdxE[Pd^)] Cum, (Tr[lF,(l)]) ; 

in) Cumfc (A/jpj) = Ylx\-k dx Cum;(,)(P) Cum, . 

iv) FCumfc (A/^pj) = Y^xhk dx Cum;(,)(P) Cum, (Tr[lFrf(l)]) . 

Proof. Moments i) follow from (14.,ip by replacing 7 with p and a with Ud-f. Factorial moments 
ii) follow from (14.3p by replacing 7 with p and a with ujd. Cumulants Hi) follow from ()4.3I) by 
replacing 7 with x-P and a with Wrf./ 3 . Note that moments of x-P are cumulants of p. Factorial 
cumulants iv) follow from (14.3p by replacing 7 with x-P and a with ojd. Q 














To compute cumulants of A/'[p], cumulants of are necessary. They can be recovered from 
(j4.5p setting rrij = m. 


5 Photocounting statistics 


When facing with sampled photocounters, two problems need to be solved: to check if the usual 
hypothesis of semi-classical theory of statistical optics hold, that is to infer about the Poisson 
distribution of photocounting, and to estimate its intensity field. Both tasks can be performed 
by using 17-statistics [9], in a different way depending on which kind of information have been 
sampled. 

If we have a random sample of overall photocounters iVt^l = ,... for each 

wave, the first task may be performed simply by checking the additivity property in (j8.5[) and so 
by estimating cumulants. Cumulants and their products can be estimated from a random sample 
X = {xi, X2, ■ ■ ■, Xn) by using a family of [/-statistics k\{x) called polykays [9]. If {cj} is the 
sequence of cumulants of iVl^l (or A/[p]), then E[k\{x)] = .... Polykays up to order 3 are: 


Ki{x) 

ki,2{x) 


£i 

5 

n 


Si — S2 

- rp-, 

n[n — 1) 

+ {n + l)siS 2 - ns 3 
n{n — l)(n — 2 ) 


K2{x) = 
, K3{x) 


ns2- si 
n (n — 1 ) 
_ 2sf - 3 
n(n - 


, Ki3{x) = 

n{n 

siS2n — in? S3 


l)(n- 2 ) 


3siS2 -|- 2 S 3 
-l)(n- 2 ) 


where Sj = power sum symmetric polynomials in the sample x for all non-negative 

integers j. 

The single index ac’s are the fc-statistics; the multi-index k’s are the polykays. The degree is 
the sum of the subscripts, that is the integer of which A is a partition. The size n of the sample 
needs to be greater than the degree. Polykays were introduced by Fisher (see [9] and references 
therein for more details) as “inherited on the average”, a property which gives to these functions 
a common interpretation independent of the sample size m- The “inheritence” property states 
that if y is a sub-sample of x obtained by simple random sampling, then E[K\{y)\x] = kx{x). 

Like cumulants, polykays enjoy of the additivity property (|3.5I) . 


Proposition 5.1. Polykays computed on the overall photocounting effect of p waves linearize in 
polykays computed on the overall photocounting effect of a single wave. 

Also equation (13.4p may be useful as a first step to verify if the underlying stochastic model is 
of Poisson type. Indeed, if photocounters have been sampled for a fixed vector of intensities (for 
example by using Monte-Carlo methods), then single polykays conditioned to known intensities 
should result approximatively constant. Moreover, Proposition ^.lOl shows that polykays are useful 
to compute Cum^ (A/[p]) and FCum^ (A/[p]). 

As mentioned at the beginning of the section, a different task consists in estimating the inten¬ 
sity field parameters from sampled photocounters. According to Proposition 13.21 this estimation 
can be easily carried out using [/-statistics for factorial moments, [/-statistics of factorial mo¬ 
ments can be recovered by expressing factorial moments {fk} in terms of power sums Sj of the 
sampled overall photocounters Again the symbolic method of moments helps in finding their 
expression, as the following proposition shows. 
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Proposition 5.2. The U-statistic fk for the k-th factorial moment o/A/Jp] is 

fk = ^'^dxsi(x)C\ with CA = - 1)!]'’*. 

X\-k 

Proof. From (I4.3p . factorial moments of a are E[{a.x)^] = X^Ahfc ^Aaz(A)CA with cx the sequence 
of cumulants of x, associated to the partition A. The result follows since can be estimated by 
sample moments Sk/n and cumulants of the singleton umbra are {(—— 1)!}, see [12]. □ 

A completely different scenario arises when we wish to predict photocounting effect from 
sampled intensities by using equation (13.ip . In this case, we need to estimate cumulants of 
Wishart random matrices. The reasons are twofold. Factorial cumulants of overall photocounters 
are equal to cumulants of Wishart random matrices and linearize on outers product of {Aij}, 
allowing to check the independence property. Moreover, since complete Bell polynomials are easily 
recovered from any symbolic packages, these estimators allow us to recover also an approximation 
of probability distribution (13.111 . Let us underline that in the literature on photocounting, factorial 
cumulants are calculated from moments of the recorded photon counts by using the classical 
moment conversion equations. Here we propose a different strategy relied on spectral polykays 
kx (A), introduced in [10]. Spectral polykays are unbiased estimators of cumulants of random 
matrices A. For Wishart random matrices, spectral polykays still have the inheritance property 
and estimate cumulants normalized to the dimension [5]: 

E{kx[Wd{p)\} = -^Cumx[Tv {Wd{p))]. (5.1) 

In difference from polykays, referring to random samples of a population, spectral polykays involve 
spectral samples. A spectral sample is the eigenvalue vector e of a random matrix A, if the size 
of the sampling n is equal to the order of the matrix A. If n < d a suitable subsample of spectral 
decomposition is selected. In the following, we assume n = d. 

Spectral polykays up to order 3 are: 

Tr[lFd(p)] _ ,, d(Tr[W,(p)])2-Tr[Wd(p)2] 

— d —’ - - dW^) - 

d(TV [Wd{p)])^-Tr[Wd{p)^] 
d(d2 - 

(Tr [Wd{p)]f{d^ - 2) - 3dTV [Wd{p)] TV [Wd{p?] + 4Tr [Wd{pY] 

d {d? - 1) {d? - 4) 

-2 dTV [Wd{pf] + {d^ + 2) Tr [Wd{p)] TV [Wd{pf] - d (TV [Wd{p)]f 

d [d? - 1) {d? - 4) 

^ 2(Tr [Wd{p)]f - 3dTr [Wd{p)] TV [Wd{pf] + d^Tr [Wd{pf] 
d{d?-l){d?-A) 

For completeness, the general formula to recover spectral polykays in terms of intensities is pro¬ 
vided in Theorem 15.31 The proof is given in [T^ . The algorithm implementing this formula is 
available in [7]. Permutation^^ of cycle structure A are involved. 

®A permutation cr of [fc] can be decomposed into disjoint cycles C{g). The length of the cycle c £ C{g) is its 
cardinality, denoted by [(c). The number of cycles of cr is denoted by |C(cr)|. Recall that a permutation cr with ri 
1-cycles, 2-cycles and so on is said to be of cycle class A = (1’’^, 2’"^, ■ ■ ■) h fe. 


(5.3) 


Ki(e) = 
K2(e) = 

ki3{e) = 

«i,2(e) = 

%(e) = 


10 










Theorem 5.3. IfXhk, then E {kx [Wd{p)]} = UjU'-Y^ Erc.=aMld)-YT)E{%x[Wd{p)] (w)} , 
where Tr(A)(cT) = ncGC(o-) with A either the matrix identity Id either the Wishart ran¬ 

dom matrix Wd{p)i and Tr(74)“^ is the inverse function^ of%x{A). 


6 Multivariate photocounting effect 

A detailed description of photocounting involves the computation of joint moments and joint 
cumulants 

(<.tl ■ ■ ■ = ™L'' Cum, ,,,, (6.1) 

with k = {ki,..., kd) € Nq and = (iVijp],..., A'djp])- In order to deal with sequences (16.ip . 
the symbolic method of moments has been generalized to multi-index k and vectors of umbral 
monomials (8 . Vectors of umbral monomials correspond to correlated random vectors when 
their supports 3 are not disjoint. Following the notations introduced in [8], denotes the 

moment generating functiorU of N'^''. In the following, for brevity, we referred to the number p of 
superimposed waves only when necessary. 


6.1 Multivariate moment symbolic method 


Let {vi,... ,Vd} € C[yl] a set of umbral monomials with support not necessarily disjoint. A com¬ 
plex sequence {ofc}, with = aki...kd = 1 , is represented by the d-tuple u = (i^i,..., Ud) 

iff 


E[u^] = Ofc, fee ni 


( 6 . 2 ) 


If {ui,... , 1 'd} are umbral monomials with disjoint supports then The 

elements in (j 6 . 2 p are called multivariate moments of v and, by analogy with random vectors, 






i>l \k\=i 


is the moment generating function of u, with z = (zi,..., Zd), |fc| = ki + - ■ -Ekd and k\ = kil • • • kdl 
Two umbral d-tuples ui and 1^2 are said to be uncorrelated if and only if E[izf i/^] = E[i/f]E[i/ 2 ] 
for all k,j G Ng. They are said to be similar if E[i/f] = E[iz 2 ] for all k G Ng, in symbols Ui = 1 / 2 . 
As done for the univariate case, if the sequence {ofc} is umbrally represented by u, its sequence 
{cfc} of multivariate cumulants satisfies 

E E z)-l]. (6.3) 

i>l \k\=i 


Next definition generalizes Definition 14.11 to multi-index partitions. 

®The inverse function f~^{a) of f{a) satisfies /b) /”^('^) = u:=a f~^ fi^) = ^('^) where &{a) = 1 

if a is the permutation identity, 0 otherwise. 

^The support of an umbral polynomial p € R\A\ is the set of all umbrae which occur. 

®Note that in the literature, the probability generating function of is erroneously called moment generating 
function since factorial moments are usually recovered from it. 
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Definition 6.1. If A is a multi-index partition and {ofe} (respectively {cfc}) is the sequence 
of multivariate moments (respectively multivariate cumulants) of v, the product a\ = ■ ■ ■ 

(respectively c\ = ‘ ‘) is said associated to the partition A. 

For example a multi-index partition A of fc = (2,1, 5) is A = (Ai, A 2 , A 3 ) with A'^ = (0,0,1), A 2 = 
( 1 , 0 , 1 ), Ag = ( 1 , 1 , 2 ) and a\ = aooi oio 1 12 - If in the dot-product 7 .a, the umbra a is replaced 

by the d-tuple u, then equation ()4.3I1 generalizes in 

^ gi(^x) cx, (6.4) 

where cx is the product of multivariate cumulants of u associated to A. 

Definition 6.2. The r.v. with sequence of moments (16.41) is named multivariate generalized 
random sum. 

Suitable choices of 7 and iz correspond to suitable choices of {g^} and {cfc} in ()6.4p and allow 
us to recover moments of special (auxiliary) umbrae. A special dot-product is the i/-cumulant 
umbra representing the sequence of cumulants {cfc} in (|6.4I) . The following definition states 
when a d-tuple of umbral monomials may be replaced by a random vector. 

Definition 6.3. A d-tuple u of umbral monomials represents a random vector X, if iz umbrally 
represents the sequence of multivariate moments {E[X^]), that is E[i/^] = E[X^] for all k G Nq. 

As in the univariate case, the dot-product y.iz represents a sum of random vectors indexed 
by a not necessarily univariate integer-value r.v, what we have called multivariate generalized 
random sum. In particular 7 . 1 / represents the iz-cumulant umbra. 

6.2 Computations of joint photocounters 

Complete Bell polynomials Yk{xi.,... ,Xk) in ()3.1I) are 

k 

Yfc(xi,... ,Xfc) = ^ ^ dxx[^xl^--- (6.5) 

i=l X\-k,l{X)=i 

Joint moments and joint cumulants can be computed by using suitable generalizations 

of complete Bell polynomials Yk{xi,... ,Xk) with the indeterminates {xi,... ,Xk} replaced by 
umbrae. More precisely, let us consider the auxiliary umbra 7 ./?.a, that is the summation 7 times 
of 13.a. The auxiliary umbra 13.a represents a compound Poisson r.v. of parameter 1, that is a 
summation N times of a r.v. represented by a, with N ~ Po(l). Moments of y./J.a, computed 
by means of equation (14.31) . result to be a first generalization of complete Bell polynomials (16.51) . 
Indeed, when a is replaced by f3.a in equation (|4.3p . we have [9] 

E[(7./3.a)^] = '^dx ggx) a\ ( 6 . 6 ) 

Ahfc 

® A partition of a multi-index A h m is a matrix A = [Xjt) of non-negative integers and with no zero columns 
in lexicographic order such that Aji + Xj 2 -!-•••= ruj for j = 1,2,... ,d. As for integer partitions, the notation 
A = (Aj^^, Aj^,...) means that in the matrix A there are n columns equal to Ai, r 2 columns equal to A 2 and so on, 
with Ai < A 2 < • ■ •. The multiplicity of Ai is Vi and we set m(A) = (ri,r 2 ,...). The number of columns of A is 
denoted by ((A). 
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with {uk} and {gk} umbrally represented by the umbra a and 7 respectively. The generalization to 
the multivariate case of equation (|6.6I) is obtained by replacing a with the d-tuple u = {vi,, I'd) 

E[(7./3.*z)^] = g ax, (6.7) 

where a\ is the product of multivariate moments of iz associated to A. More general expressions 
of equation (16.7p correspond to moments of with {izi,..., Vd} d-tuples of 

umbral monomials [ 8 ]. Set 

(71./3.1/1 -I-h jd-P-i'd)’^ = ■ ■ • )7d)- (6.8) 

The polynomials ’*^'^^ 71 , • • • j 7 d) are said generalized complete Bell polynomials. For d = 1 

we recover yj^\'y) = {y.fi.u)^. By using a suitable generalization of multinomial expansion [ 8 ], 
the fc-th moment of ( 71 ./3. 1/1 7 d./ 3 .iZrf) is 


E 




(7i,---,7d) 


E 




k 

; ■ ■ ■ ; '^d 


E 


(7i./3.izi)*i • • • (7rf./3.izd)*'^ . (6.9) 


In the following, let us denote by ij the umbra representing the intensity Ij for j = 1,... ,d and 
by L the corresponding d-tuple l = {li, ... ,id). 


Theorem 6.4. If k £ Nq, then mk = E ^ t^) 


with 


Uj = {£,... ,. ..,£), for j = 1 ,... ,d. 

J-th place 

Proof. From (I6.10p . we have f{Uj,z) = for j = 1,2,... ,d and 


( 6 . 10 ) 


. i=^ 


[=E, 

f \ 

exp 




^t,.{/(u„2) - 1} 

i=i 


( 6 . 11 ) 


The result follows by observing that the right hand side of (|6.11l) is the moment generating 
function of ii .fl.Ui + ■ ■ ■ + id -fl-Ud. □ 


Next corollary gives the explicit expression of joint moments of photocounters. 

Corollary 6.5. If k = (ki,..., kd) € Ng then mk = E {n7i [EliiSft.ij'j]}. 

Proof. Let us consider the expansion (|6.9p with (iz^,..., Ud) replaced by (wi,..., Ud) and ( 71 ,..., 7 ^) 
replaced by (ti,..., t^). Since 


E 


ui 

=j 

i 

1 


0, if j is such that jfc / 0 for A; = 1, 2,... , i — 1, z -|- 1,... d, 
1 , otherwise 


( 6 . 12 ) 


among the vectors (ii,... ,id) giving not zero contributions in (j6.9p . there are those satisfying 
(ij). = kjhj^i for j = 1, 2,... , d. In this case (^^ ^ =1. Since 


E 


(ii ./3.iii)*i • • • (td =e|e (ii ./3.rii)*i • • • (zd ./3.Ud)*'^ I ii,... ,Zd | 


(6.13) 


the result follows by multiplying E[(z,j ./3.Uj)7 | ij] = E[(ij ./3)^z = Ylt=i i ~ 

1,2,... ,d and getting the overall expectation of the product. □ 
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Corollary 6.6. For k = {ki,..., k^) € Ng, if there exists j G {1, 2,..., d} sueh that kj = 0 then 

ruk = 0 . 


Remark 6.7. To compute mk by using Theorem 16.41 a symbolic procedure is available on de¬ 
mand. This consists in expanding the product in the right hand side of (|6.9I1 and then replacing 


occurrences of with E 


jk\ jkd 


. The algorithm nCWishart [ 6 ] allows us to compute 


E 


jki 


. Indeed, this algorithm allows us to compute general joint moments 


E 


Joint moments E 


{Tr [Wd{p) • • • Tr [W^ip) } with ,..., € 

can be recovered from (j6.14p by choosing 


'^dxd 


(6.14) 


jk^ ... 


iHi)st = 


1 , if i = s = t, 
0 , otherwise, 


s,t = 1,2,d, and i = 1,2,..., d. 


Denote by the A:-th multivariate factorial moment of N. As the following theorem states, 
also multivariate factorial moments can be expressed via generalized complete Bell polynomials. 


Theorem 6.8. If k ^ Ng, then = E (t^^ ■ ■ ■ ,id) 


with 


Xj = ••,£), for j = l,...,d. 

J‘th place 

Proof. By using the same arguments employed in the proof of Theorem l6.4[ the moment generating 
function of the factorial moments of N is 


. 1=1 


[=e| 

( \ 
exp 


V 


'^i'AfiXi,w) - 1 } 

1=1 


(6.15) 


where /(Xi,re) = 1 + Wi and w = {wi,... ,Wd). The result follows by observing that the right 
hand side of (I6.15P is the moment generating function of ii .(3-Xi + ' " + I'd •AXd- D 

Corollary 6.9. If k = {ki,..., kd) G Ng, then ^k = E (l^^ ' " k 


Proof. From Theorem 16.81 and equation (j 6 . 8 p we have ffc = E {ii.ui+ ■■■ +td.Ud) 
AXj = foi j = 1,... ,d. Then by using the multinomial expansion (j6.9p . we have 


E 


{il.Ui H-h Ld-UdY 


E 


k 
? • • • ; 


E 


(il.Ul)*l • • ■ (Ld-UdY'^ 


Since 


(6.16) 


In order to evaluate products on the right hand side of (I6.16[ 


Since ij-Uj = ij.j3.x-Uj and = Xji in evaluating E 

A ■■ ■ YY gives contribution, from which the result follows. 

Joint cumulants can be computed by using the additivity property on waves. 


, equation (|6.7p has to be employed. 
nj=i(b-Wj)*i only joint products 

□ 
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Theorem 6.10. If k £ Nq, then 


p 

Cumfc(ArW) = k\J 2 Y^ 

i=\ XGPk 


- 1 ]! 

m(A)!A! 


n® 


d 

i^s)j 


n 

^ S{{Xs)j,t)\Xij +miP^ 


_i=i 

t=i 



(6.17) 


where Pi = {X = {X'^.X "'^,...) h fc : {Xs)j 7 ^ 0, Vj = 1,..., d, s = 1, 2,...}. 
Proof. First, let us prove the additivity property 


p 

Cumfc(A/'W) = y~]Cumfc(Ar'^^), (6.18) 

i=l 


where Cumfc(Arf^') is the A:-th cumulant of the multivariate photocounter of the i-th wave with 
i = Indeed, denote by 7 *^ the umbral monomial representing \Xij + in (j2.ip . 

Then ti.fl.Ui + • • • + = X)j=i(7ii + ''' + 'ypj)-P-Uj. Since for fixed j the umbral mono¬ 

mials {71 j,..., 7 pj} are uncorrelated, equation (|6.18p follows by observing that X)j=i( 7 ii + 

- 1'rpj)-f3-Uj = H- ^Ipj-P-Uj) = Ef=i(7ii-^-Wi H- ^-fid-fl.Ud). Since 

{'jii./I.ui 'jid-P-Ud) denotes the multivariate photocounter from d pixels hit by the i-ih. 

wave with i = 1 ,... ,p, the result follows from the additivity property of multivariate cumulants. 
To get equation (I6.17h . the explicit expression of Cumfc(Ali^^') has to be computed by evaluating 
the fc-th moment of 7.(74 i./3.rii ^id-fl-Ud). The fc-th moment of 7 . 1 / is [ 8 ] 


nM’^] 


^ k\ 
^m(A)!A! 


(-l)'W-n/(A)-l]!aA, 


with a\ the product of multivariate moments of u associated to A. The result follows by replacing 
a\ with the product of multivariate moments of 'jid-P-Ud) associated to A. 

Moments of {■jn.fl.ui 'yid-(3-Ud) are given in Corollary 16.51 with tj replaced by 7 *j. □ 

Within estimation, the importance of Theorem 16.101 relies on the circumstance that if estima¬ 
tors of Cumfc(ArM) linearizes, according to (|6.18p the underlying distribution of photocounters 
can be assumed of mixed multivariate Poisson type. Unbiased estimators of Cumfc(Ar)^') and 
Cumfc(Ar[^') can be computed using multivariate polykays [9]. 

Theorem 6.11. The multivariate Poisson-Mandel transform admits the following expansion in 
series 


p(Ar = k) = 


kl 


, i>i \j\=i j • V i>i \j\=i 


(l)b1 


E 


U + ky. 


j\ ^ tn(A)!A! 

■’ Ahfe+i ^ ' 


c\ 


with c\ the product of multivariate cumulants of I associated to A. 
Proof. The first equality follows from (| 1 . 1 |) by observing that 


exp{-(/i -h • • • -F 7rf)} = 1 -F ^ ^ 

*>i \j\=i 


,.(-1)1^1 


j! 


The second equality follows from equation (16.41) with gpx) replaced by 1, since u = (/3.7).iz. □ 
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If sampled intensities are available, multivariate polykays for I can be employed to approximate 
c\ in Theorem 16.111 and also factorial cumulants of Af, as the following theorem shows. 

Theorem 6.12. FCumfe (AT) = Cfc. 

Proof. Denote by K{N, z) the cumulant generating function of Af, with K(N, z) = log[f{N, z)]. 
The factorial cumulant generating function is 

K[N, z) I Zi = log(l + u;P = log[/ (Af, z)] I Zi=log(l + u;P = log[/ (t, Z)], 

i=l,2,...,d i = l,2,...,d 

where the last equality follows from Corollary 16.91 □ 

7 Conclusions and open problems 

This paper has introduced the symbolic method of moments as an efficient tool to deal with photon 
statistics. Instead of using factorial moments, as usually proposed in the literature, cumulants 
have been employed due to their additivity property which simplifies the inference on the process 
as well as factorial cumulants. Moreover these sequences have a simpler expression for Wishart 
random matrices whose traces are the intensities of photocounting. In addition, the algorithms 
available for computing unbiased estimators of cumulants, that are polykays and multivariate 
polykays, have efficient implementation. 

Several open problems arise from the methodology suggested in this article. As investigated 
in m, there is a connection between Wishart random matrices and natural exponential families. 
So distributions of photocounters could be studied by using properties of natural exponential 
families. A first attempt in this direction has been presented in [3] and [l3], where only special 
cases are considered. A more general development, involving cumulants, could turn to be useful in 
finding sufficient photon statistics relied on maximum likelihood method. The use of the symbolic 
method of moments within these applications is currently under investigation. 

A further development consists in studying photocounting vs. time, through a compound 
Poisson stochastic process. The so-called compensated version consists in normalizing the process 
with respect to its temporal mean. Cumulants play a fundamental role in dealing with compen¬ 
sated Poisson process [23] and we believe that the techniques here introduced should be fruitfully 
applied. 
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